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SUMMARY 

The  cohesive  element  approach  is  getting  increasingly  popular  for  simulations  in  which  a  large  amount 
of  cracking  occurs.  Naturally,  a  robust  representation  of  fragmentation  mechanics  is  contingent  to  an 
accurate  description  of  dissipative  mechanisms  in  form  of  cracking  and  branching.  A  number  of  cohesive 
law  models  have  been  proposed  over  the  years  and  these  can  be  divided  into  two  categories:  cohesive  laws 
that  are  initially  rigid  and  cohesive  laws  that  have  an  initial  elastic  slope.  This  paper  focuses  on  the  initially 
rigid  cohesive  law,  which  is  shown  to  successfully  capture  crack  branching  mechanisms  in  simulations. 
The  paper  addresses  the  issue  of  energy  convergence  of  the  finite-element  solution  for  high-loading  rate 
fragmentation  problems,  within  the  context  of  small  strain  linear  elasticity.  These  results  are  obtained  in 
an  idealized  one-dimensional  setting,  and  they  provide  new  insight  for  determining  proper  cohesive  zone 
spacing  as  function  of  loading  rate.  The  findings  provide  a  useful  roadmap  for  choosing  mesh  sizes  and 
mesh  size  distributions  in  two  and  three-dimensional  fragmentation  problems.  Remarkably,  introducing  a 
slight  degree  of  mesh  randomness  is  shown  to  improve  by  up  to  two  orders  of  magnitude  the  convergence 
of  the  fragmentation  problem.  Copyright  ©  2006  John  Wiley  &  Sons,  Ltd. 
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1.  INTRODUCTION 

Structures  often  develop  complex  fracture  and  fragmentation  patterns  during  the  failure  process, 
and  the  quest  for  understanding  the  pattern  and  distribution  of  fragment  size  has  captured  the 
attention  and  imagination  of  the  scientific  community  for  many  decades.  Models  for  predicting 
fragment  size  distribution  have  been  applied  to  understand  the  physics  of  hypervelocity  impact 
[1,2],  crash  performance,  explosive  drilling  [3],  and  even  the  size  distribution  and  clustering  of 
galaxies  resulting  from  the  formation  of  the  universe  during  the  big  bang  [4].  Armor  ceramics 
are  another  important  class  of  materials  that  fail  on  impact  due  to  a  fast  fragmentation  pro¬ 
cess  [5].  It  has  been  observed  that  the  interaction  between  the  penetrator  and  the  fragmented 
zone  ahead  of  it  (often  referred  to  as  the  comminuted  or  Mescal  zone)  is  key  to  penetration 
performance  [6]. 

Despite  this  large  body  of  work,  our  understanding  of  the  physics  of  fragmentation  is  still 
incomplete.  The  difficulty  emanates  from  the  inherent  multi-scale  nature  of  the  fracture  process. 
In  the  case  of  ceramic  materials,  it  is  well  known  that  strength,  crack  initiation,  and  crack  prop¬ 
agation  are  affected  by  the  presence  of  flaws  at  the  micrometer  scale.  Under  dynamic  loading 
conditions,  cracks  will  initiate  at  these  flaws,  and  potentially  propagate  catastrophically  to  cause 
large-scale  structural  failure.  Multiple  cracks  will  initiate  at  seemingly  random  locations  and  ma¬ 
terial  failure  will  occur  through  a  complex  communication  process  of  stress-wave  interactions 
between  cracks. 

Evidently,  due  to  the  space  and  time- varying  nature  of  the  problem,  the  use  of  classical  fracture 
mechanics  analysis  is  limited.  Recent  developments  in  the  field  of  computational  solid  mechan¬ 
ics,  including  large-scale  parallel  simulations,  enable  a  fresh  look  at  fragmentation  mechanics. 
Novel  simulation  techniques  permit  a  new  insight  and  an  ever-closer  agreement  between  ex¬ 
periments  and  modelling.  An  example  of  such  a  numerical  approach  is  the  cohesive  element 
method.  Cohesive  elements  originate  from  the  concept  of  cohesive  zone,  which  was  first  in¬ 
troduced  by  Dugdale  [7]  and  Barrenblatt  [8].  The  implementation  of  a  cohesive  zone  into  nu¬ 
merical  analysis  takes  the  form  of  cohesive  elements,  which  explicitly  simulate  the  crack  pro¬ 
cess  zone.  Many  examples  of  applications  that  utilize  cohesive  elements  may  be  found  in  the 
literature  [9-22]. 

The  cohesive  element  approach  has  gained  in  popularity  over  the  years.  This  is  principally 
due  to  its  ease  of  implementation  and  the  clear  physical  picture  that  is  given  by  the  explicit 
representation  of  cracks.  Crack  branching  and  fragmentation  are  often  presented  in  the  literature 
as  natural  outcomes  of  the  method.  Nonetheless,  despite  the  success  of  the  approach  several  topics 
of  controversy  remain.  A  first  item  that  will  be  briefly  discussed  in  the  paper  is  the  choice  of  the 
cohesive  law.  We  have  shown  in  a  different  work  that  the  shape  of  the  cohesive-law  softening 
curve  has  little  effect  on  the  fragmentation  results  [20,  23].  This  paper  illustrates  that  initially  rigid 
cohesive  laws  can  successfully  model  the  ‘side-branching’  fracture  mechanism  [24,25].  Perhaps 
a  more  fundamental  issue  has  to  do  with  convergence  of  the  method.  Recently  Papoulia  et  al. 
have  provided  computational  evidence  of  mesh  convergence  of  the  crack  path.  Convergence  was 
observed  for  a  special  category  of  meshes  referred  to  as  pin  wheel-based  meshes  [26].  Crack  path 
convergence  will  not  be  addressed  in  this  paper,  which  will  be  solely  and  fully  devoted  to  the 
evolution  of  microcrack  density  as  function  of  mesh  size.  In  the  creation  of  new  surfaces,  each 
opening  cohesive  element  dissipates  a  given  amount  of  energy  (e.g.  cohesive  energy).  The  total 
energy  dissipated  in  the  cracking  problem  is  therefore  clearly  related  to  the  crack  path  and  thus 
mesh  convergence.  Energy  convergence  of  the  finite-element  solution  has  been  mainly  eluded 
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in  the  literature.  It  has  been  generally  observed  that  for  a  fixed  strain  rate,  finer  meshes  lead 
to  higher  energy  dissipation  due  to  more  microcracking  [17].  To  the  best  of  our  knowledge,  no 
demonstration  of  solution  convergence  has  been  shown  in  large-scale  fragmentation  problems. 
The  issue  of  energy  convergence  during  the  cracking  process  is  of  fundamental  importance  as  it 
clearly  affects  the  physical  robustness  of  the  numerical  methodology;  it  is  the  main  purpose  of 
the  paper. 

Section  2  of  the  paper  briefly  describes  the  finite  element  methodology,  including  the  cohesive 
element  approach.  It  also  provides  background  regarding  key  length  scales  and  time  scales  associ¬ 
ated  with  cohesive  elements.  A  two-dimensional  crack-branching  simulation  is  given  to  illustrate 
the  approach.  Section  3  analyses  convergence  issues  in  a  simplified  one-dimensional  setting  and 
implications  for  the  extension  to  multidimensional  finite  element  simulations. 


2.  FINITE  ELEMENT  METHODOLOGY,  ASSOCIATED  LENGTH  SCALES 

AND  TIME  SCALES 

We  begin  with  a  brief  review  of  the  adopted  finite  element  framework.  The  context  is  explicit 
dynamics,  small  strain  linear  elasticity.  More  details  may  be  found  in  References  [12, 14,  21],  or 
in  standard  finite  element  books  [27]. 

2.7.  Explicit  dynamic  finite -element  analysis 

We  discretize  a  body  Bq  with  line- segment  elements  in  ID,  quadratic  triangular  elements  in  2D 
(T6),  or  quadratic  tetrahedral  elements  in  3D  (T10).  Upon  discretization,  the  principle  of  virtual 
work  applied  to  the  equations  of  equilibrium  renders: 

M(p  +  Rint((p)  =  Rext(/)  (1) 

where  Rext  and  Rint  are  the  external  and  internal  forces  arrays,  M  is  the  mass  matrix,  and  <p  is  the 
nodal  co-ordinates  array. 

This  equation  is  integrated  along  the  time  axis  by  the  second-order  accurate  explicit  scheme. 
The  explicit  version  of  the  Newmark  scheme  is  obtained  taking  the  Newmark’s  parameters  to  be 
ji  =  0  and  y  =  J>,  as 

(Pn+l  =(Pn+  Atvn  +  \ At2a„ 

a;!+ 1  =  M 1  (R^j  —  R^])  (2) 

v„+ 1  =  v„  +  jAt(an  +  a„+i) 

where  v  and  a  are  the  material  velocity  and  acceleration  fields.  In  this  scheme,  the  mass  matrix 
M  is  lumped  and  diagonal. 

To  guarantee  the  stability  of  the  time  integration,  the  time  step  At  must  be  lower  than  a  critical 
value,  A  Stable,  which  is  related  to  the  dilatational  (the  fastest)  wave  speed  and  the  (smallest)  mesh 
size.  In  our  simulation,  the  time  step  is  taken  as 

At  =  CA/Stabie  —  C  min  |  J  (3) 

mesh\ya  +  2A7p/ 
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where  C  is  a  time  step  scale  factor  (C<1.0,  an  ordinary  value  is  about  0.1),  hQ  is  the  dimension 
of  the  element,  A,  p,  and  p  are  the  Lame  coefficients,  and  the  density  of  the  material.  The  quantity 
y/(A  +  2 p)/p  is  the  longitudinal  wave  speed  in  an  unbounded  medium.  In  a  ID  problem  the 
interpretation  of  hQ  is  straightforward;  in  multidimensional  setting  the  radius  of  the  in  circle  or  in 
sphere  are  taken  as  the  element  sizes. 

2.2.  Cohesive  element  method 

In  the  finite  element  analysis,  the  specimen  is  discretized  with  ordinary  continuum  elements. 
The  interfaces  between  two  neighbouring  elements  (line  segments  in  2D  or  facets  in  3D)  are 
treated  as  possible  sites  for  cracks;  see  Figure  1(a)  for  a  2D  schematic.  We  use  the  initially 
rigid  cohesive  law,  Figure  1(b),  proposed  by  Camacho  and  Ortiz  [12]  and  at  a  later  time  by 
Pandolfi  et  al.  [14]  in  a  3D  setting.  Using  this  type  of  law  implies  that  cohesive  elements  have 
to  be  added  on  the  fly,  a  scheme  sometimes  referred  to  as  dynamic  insertion  [28].  The  initial 
finite-element  mesh  is  free  of  cohesive  elements,  and  as  the  dynamic  simulation  proceeds,  co¬ 
hesive  elements  are  inserted  at  locations  where  the  stress  exceeds  a  critical  value  oc.  Several 
methods  have  been  proposed  for  identifying  what  constitutes  a  sensible  insertion  criterion.  A  dis¬ 
cussion  on  various  methods  for  dynamically  inserting  cohesive  elements  into  the  finite  element 
mesh  is  given  by  Papoulia  et  al.  [29].  The  two-dimensional  results  shown  in  this  paper  were 
obtained  by  using  an  average  of  the  normal  component  of  the  stress  tensor  to  an  element  edge. 
For  quadratic  elements,  four  Gauss  Points  contribute  to  this  average.  Clearly,  the  implementa¬ 
tion  of  a  dynamic  insertion  scheme  necessitates  extensive  computational  bookkeeping,  but  this 
methodology  avoids  artificially  altering  the  elastic  properties  of  the  medium  or  structure.  Indeed, 
in  other  cohesive  law  models,  such  as  the  exponential  cohesive  law  [24],  the  element’s  initial 
response  amounts  to  adding  a  spring  in  the  mesh,  and  this  may  substantially  change  the  problem’s 
response. 

The  initially  rigid  cohesive  model,  which  is  used  in  the  remainder  of  the  paper,  is  illustrated 
in  Figure  1(b).  We  have  shown  in  previous  work  [23]  that  the  shape  of  the  unloading  curve  does 


Bulk  or  vol umatric  element 
Cohesive  element 


Figure  1.  (a)  Finite  element  discretization  with  triangular  bulk/continuum  elements  (six  nodes)  and 
corresponding  cohesive  elements;  and  (b)  initially  rigid,  linear-decaying  irreversible  cohesive  law. 
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not  change  significantly  the  fragmentation  results.  Therefore,  we  use  a  simple  irreversible  linear 
decaying  cohesive  law: 


—  =  1  —  —  for  (5>0,  d  =  3max 
&C  3C 


t  _  3 

G  c  ^max 


—  for  a  <  a  max 

3c 


(4) 


where  t  is  an  effective  traction,  3  is  an  effective  displacement,  oc  is  the  maximum  cohesive  force, 
and  3C  is  the  critical  opening  distance.  Note  that  a  cohesive  element  resists  opening  until  it  is 
fully  damaged  (<5max  =  t)c).  The  notion  of  an  effective  displacement  should  be  further  explicated. 
The  original  idea  developed  by  Camacho  and  Ortiz  [12]  is  that  the  response  of  a  cohesive  element 
should  depend  on  a  combination  of  normal  and  shear  deformations.  They  defined  the  effective 
opening  displacement  3  as 


(5) 


where  <5n  and  3S  are  the  normal  and  shear  opening  displacements  over  the  cohesive  surface.  The 
parameter  jS  assigns  different  weights  to  the  sliding  and  normal  opening  displacements.  Similarly, 
the  effective  traction  t  is  defined  as 


t 


(6) 


where  ts  and  tn  are  the  shear  and  normal  tractions.  Clearly,  a  negative  3n  indicates  that  contact  is 
occurring.  In  this  paper  we  enforce  a  simplistic  contact  algorithm  in  which  contacting  facets  are 
modelled  as  a  perfect  junction  with  no  overlap.  At  contact  locations,  compressive  stress  waves 
travel  through  the  material  as  if  it  were  undamaged  until  tensile  waves  reopen  these  cohesive 
elements. 

Three  material  parameters  are  associated  with  this  choice  of  cohesive  law.  The  physical  meanings 
of  oc  and  3C  have  been  described  earlier.  The  third  parameter,  which  is  not  independent  of  the 
other  two  but  is  often  easier  to  evaluate,  corresponds  to  the  area  under  the  curve  of  Figure  1(b). 
It  is  the  fracture  energy  that  is  needed  to  fully  open  a  unit  area  of  crack  surface: 


Gc  =  2rc  = 


oc3c 

~Y~ 


(7) 


where  Yc  is  the  surface  energy. 

Clearly,  the  success  of  simulations  using  cohesive  elements  resides  largely  in  the  proper  choice 
of  the  material  parameters.  Additionally,  it  is  contingent  to  resolving  appropriate  length  scales  and 
time  scales  associated  with  the  fragmentation  process.  The  next  section  revisits  some  key  concepts 
in  this  regard. 


2.3.  Associated  length  scales  and  time  scales 

Various  length  scales  and  time  scales  are  associated  with  cohesive  elements.  They  are  functions 
of  elastic  and  fracture  material  properties.  The  following  convention  is  adopted:  E  is  the  Young’s 
modulus,  p  the  mass  density,  c  the  wave  speed,  oc  the  critical  stress  at  which  the  decohesion 
process  begins,  Gc  the  critical  energy  release  rate,  3C  the  critical  opening  displacement,  and  s  the 
loading  rate. 
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2.3.1.  Cohesive  zone  concept.  In  order  to  properly  capture  the  physics  of  dynamically  failing 
material  in  the  vicinity  of  the  crack  tip  in  numerical  simulations,  it  is  important  that  the  finite 
element  mesh  size  be  smaller  than  the  cohesive  zone  length.  The  number  of  finite  elements  needed 
to  model  the  behaviour  at  the  tip  of  the  crack  will  be  highly  dependent  upon  the  constitutive 
description  of  the  cohesive  zone  and  bulk  continuum,  and  the  nature  and  rapidity  of  the  internal 
and  external  loadings. 

If  the  cohesive  zone  tractions  or  ‘failure  stresses’  are  constant  and  independent  of  position  along 
the  failure  zone,  viz,  o  f(£)  =  oc,  where  £  is  the  co-ordinate  along  the  cohesive  zone,  a  number  of 
authors  [30-32]  have  shown  that  the  cohesive  zone  length  a  can  be  written  as 


7i  / K\\2  7i  EGC 
8  \  oc  )  8  cr2(l  —  v2) 


(8) 


where  K\  is  the  mode-I  stress  intensity  factor.  The  expression  for  the  cohesive  zone  length  is 
obtained  by  solving  a  boundary  value  problem  using  the  method  of  superposition;  here,  it  is 
assumed  that  it  is  permissible  to  sum  the  stresses  0®  in  the  neighbourhood  of  the  crack  tip  without 
a  cohesive  zone  due  to  loads  acting  in  and  on  the  elastic  continuum,  with  the  stresses  induced  in 

f 

the  solid  due  to  a  failure  stress  distribution  oy  at  the  tip  of  a  crack  with  a  cohesive  zone  [8].  Using 
the  notation  in  Reference  [31],  the  superposition  can  be  expressed  as 


Gy  =  <7°  +  Gy  (9) 

where,  for  example,  the  y-component  of  stress  in  a  linear  elastic  medium  on  y  =  0,  and  close  to 
the  crack  tip,  is  given  by  Williams  [33]  as 


Kimo 

V^V 


(10) 


where  H  is  the  Heaviside  function,  and  the  stresses  due  to  the  action  of  the  cohesive  zone  failure 
stresses  [8, 31]  very  close  to  the  crack  tip  are  given  by 


z  =  r GfOtm 
y  Tly/Z] [  JO  VI 


+  HOT 


(ID 


Substitution  of  Equations  (10)  and  (11)  into  Equation  (9),  yields  the  resultant  stress  in  the  linear 
elastic  continuum. 


„,=  **«■>_  >  f^  +  Hor 

V^VJi  nVIl  Jo  VI 


(12) 


Since  the  higher-order  terms  in  Equation  (12)  are  finite,  then  oy  will  be  finite  if  and  only  if  the 
first  two  terms  can  be  equated,  viz, 


Ki  = 


Gf(0  d£ 

VI 


(13) 


By  substituting  the  assumption  that  g/(c)  =  gc  into  Equation  (13)  the  result  given  in  Equation 
(8)  is  readily  obtained.  In  another  case,  where  it  is  assumed  that  the  failure  stress  varies  linearly 
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within  the  cohesive  zone  <Jf(0  =  gc(1  —  £/a),  Palmer  and  Rice  [34]  (and  also  Rice  [32])  derive 
a  cohesive  zone  length,  which  is  |  times  the  estimate  in  Equation  (8),  i.e. 


9n/Ki\2_9n  EGC 
32\ac)  32  <t2(1  —  v2) 


(14) 


Equations  (8)  and  (14)  provide  static  estimates  for  the  cohesive  zone  length.  Clearly,  a  simulation 
that  includes  cohesive  elements  should  have  elements  smaller  than  this  dimension  to  resolve 
the  fragmentation  process.  This  procedure  is  now  commonly  adopted  in  the  literature  [24,35]. 
Nonetheless,  later  in  this  paper,  we  will  question  the  appropriateness  of  using  this  length  scale  in 
dynamic  fracture  and  fragmentation  simulations. 


2.3.2.  Time  scale  associated  with  opening  of  cohesive  element.  Once  a  fine  enough  finite-element 
discretization  is  obtained,  the  use  of  cohesive  elements  results  in  an  additional  length  scale.  This 
is  a  consequence  of  the  time  associated  with  the  opening  of  an  isolated  microcrack  [12].  Part 
of  the  analysis  that  follows  was  developed  in  Reference  [23]  and  only  the  essential  features  are 
summarized  herein. 

For  argument’s  sake,  we  consider  a  bar  made  of  a  linear  elastic  material.  The  bar  is  loaded  uni¬ 
formly  in  tension  up  to  rupture.  The  cohesive  behaviour  follows  an  initially  rigid  linear  decreasing 
cohesive  law  (Figure  1(b)). 

The  following  ordinary  differential  equation  governs  the  monotonic  unloading  behaviour  in  an 
isolated  cohesive  element  [23]: 


where  (5coh  denotes  the  opening  displacement  of  the  cohesive  element  (initially  it  is  zero).  The 
ordinary  differential  equation  (15)  can  be  simplified  by  normalizing  the  variables  similarly  to 
Drugan  [36]: 


T-  <5coh 

dcoh  —  5 

dc 


2  Gct 
pcdc 


t 

to’ 


(16) 


where  the  characteristic  time  to  [12],  to  be  understood  as  the  ‘response’  time  of  the  cohesive 
element,  is 


pcSc  =  E  T8C\EGC 
2gc  2gc\c  )  cg2c 


(17) 


and  to,  a  derivative  of  the  characteristic  time  scale,  is  the  characteristic  strain-rate  [23]: 

gc/E  2  g2  2cg2  cg\ 

60  =  ~7T  =  ~prc^c  =  e^c  =  ~eAg~c 


(18) 


Note  that  this  is  also  an  intrinsic  material  parameter,  and  determines  whether  a  particular 
externally  applied  strain  rate  is  perceived  as  a  low  or  high  strain  rate  during  this  dynamic  decohesion 
process.  The  non-dimensional  ordinary  differential  equation  and  initial  condition  are: 


<5coh  : 


d<5coh 

— E - £t 

d  t 


(19a) 
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and 


<5cohlF=o  =  0  (19b) 

This  equation  is  readily  solved  as 

4>h  =  e[exp(F)  -7-1]  (20) 

Equation  (20)  describes  the  dynamics  of  an  isolated  decohesion:  the  evolution  of  the  cohesive 
crack  opening  displacement  for  an  isolated  crack  in  an  expanding  bar.  The  crack  opening  dis¬ 
placement  is  a  non-linear,  monotonically  increasing  function  of  time,  controlled  by  the  material 
properties  and  the  external  loading  s.  Substituting  (20)  into  the  cohesive  law,  the  non-dimensional 
cohesive  stress  may  now  be  expressed  as 

acoh  =-™=l-  e[exp(0  -t-  1]  (21) 

The  critical  time  tc  at  which  the  decohering  point  completely  fails  can  be  computed  by  setting 
aC0h  =  0  in  Equation  (21): 

e[exp(7c)  -  tc  -  1]  =  1  (22) 

For  the  purpose  of  this  study,  we  are  interested  in  the  high  strain  rate  response,  in  which  case 
£^1,  and  tc  is  small.  Equation  (22)  reduces  to 

j 2  _  i 

[exp(7c)  -  tc  ~  1]  =  f  +  O  (tl)  =  -  (23) 

^  £ 

which  results  in 

/2\  1/2 

ic=  T  (24) 

\£  / 

Thus  the  normalized  time  for  complete  decohesion  is  inversely  proportional  to  the  square  root 
of  the  external  strain  rate  for  very  high  strain  rates.  In  dimensional  form  this  may  be  written  as 


This  time  scale  is  of  practical  importance.  In  order  to  properly  resolve  the  unloading  part  of  the 
cohesive  law  in  a  dynamic  fragmentation  simulation,  the  integration  time  step  must  be  roughly  an 
order  of  magnitude  smaller.  We  will  come  back  to  this  point  in  view  of  numerical  results. 

2.4.  Illustration  of  cohesive  element  approach:  crack  branching  simulations  in  PMMA 

Before  discussing  the  topic  of  convergence,  we  illustrate  the  cohesive  element  approach  with  results 
of  dynamic  crack  propagation  and  branching  in  PMMA.  Crack  branching  tests  are  increasingly 
becoming  a  benchmark  test  and  a  number  of  recent  papers  have  used  it  for  validation  purposes 
[21, 24,  35,  37,  38].  Here,  the  simulations  of  Falk  et  al.  [24]  were  reproduced  to  test  the  ability  of  the 
initially  rigid  cohesive  elements  to  capture  crack  branching  in  a  2D  setting.  Previously,  we  obtained 
probing  crack  branching  results  in  3D  PMMA  plates  [21].  Nonetheless,  it  had  been  observed  that 
the  type  of  cohesive  law  used  in  this  paper,  might  fail  to  capture  branching  mechanisms  [24]. 
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Figure  2.  (a)  Finite- element  mesh  o  a  PMMA  plate  under  mode  I  opening  (the  mesh  contains  36  096  T6 
elements);  and  (b)  simulation  results  showing  crack  branching,  which  were  obtained  for  Wo/Gc  =  13.5 
(vertical  displacement  contours  are  shown;  displacements  on  top  and  bottom  boundaries  are  ±43  pm). 


In  order  to  test  the  approach  in  2D  as  well,  we  have  developed  a  2D  dynamic  insertion 
scheme  for  cohesive  elements.  Similar  to  Falk  et  al  [24]  we  meshed  a  h  =  3  mm  square  block 
of  PMMA  with  an  edge  crack  extending  0.25  mm  into  the  block  from  the  left  edge.  The  elastic 
properties  for  PMMA  are  chosen  to  be  £  =  3.24  GPa,  v  =  0.35,  and  p  =  1.19  g/cm3,  while  the 
fracture  properties  are  Gc  =  300  N/m,  and  oc  =  130  MPa.  Finite-element  meshes  were  obtained 
for  five  different  mesh  sizes  hQ  =  64,  48,  40,  32,  and  12  pm,  corresponding  to  a  total  number 
of  nodes  of  10211,  18201,  26217,  41029,  and  72713.  Figure  2(a)  represents  the  finest  mesh 
(72713  nodes,  36096  elements).  The  cohesive  zone  length  estimate  obtained  with  Equation  (14) 
is  58.5  pm.  Since  the  cohesive  zone  size  is  bigger  than  all  but  the  coarsest  mesh  size  used  in 
this  study,  we  expect  to  be  able  to  resolve  the  important  details  of  the  dynamic  crack  propagation 
process. 

Different  loading  conditions  were  tested  including  the  ones  chosen  in  Reference  [24]  and  in 
Reference  [37].  In  both  cases,  the  loading  is  applied  uniaxially  by  prescribing  a  displacement 
in  the  vertical  direction  on  the  top  and  bottom  boundaries  (the  side  boundaries  are  stress  free). 
In  Reference  [24],  the  displacement  is  applied  dynamically,  starting  from  uniformly  applied  velocity 
gradients  in  x  and  y  directions,  and  thereafter  maintaining  a  constant  strain  rate  in  y  direction 
on  the  top  and  bottom  boundary.  This  form  of  loading  does  not  generate  waves  propagating  in 
from  the  boundaries  but  creates  a  constant  source  of  energy  input  in  the  propagating  crack  tip. 
While  branching  results  were  obtained  for  all  meshes  with  this  form  of  loading,  we  preferred  the 
boundary  conditions  detailed  in  Reference  [37].  There,  we  prescribe  a  fixed  vertical  displacement, 
<5o,  on  the  top  and  bottom  boundary  and  solve  for  the  static  problem.  The  strain  energy  per  unit 

Copyright  ©  2006  John  Wiley  &  Sons,  Ltd.  Int.  J.  Numer.  Meth.  Engng  2007;  69:484-503 

DOI:  10.1002/nme 


THE  COHESIVE  ELEMENT  APPROACH  TO  DYNAMIC  FRAGMENTATION 


493 


Figure  3.  Cohesive  energy  dependence  on  mesh  size  for  crack  branching 
simulations  in  PMMA  plate  for  Wq/Gc  =  13.5. 


length  in  a  plate  without  a  crack  under  fixed-grip  boundary  loading  is  well  defined: 


W0  = 


2ES 

~~h 


2 

0 


(26) 


Wo  becomes  available  for  creating  new  surfaces  when  the  simulation  begins  (experimentally  this 
would  be  equivalent  to  introducing  a  sharp  crack  with  a  razor  blade  in  a  pre-strained  plate  [39]). 
Naturally,  Wo  is  a  direct  function  of  the  applied  displacement  field.  The  larger  the  displacement 
the  more  energy  is  available  for  propagating  the  crack.  At  sufficiently  high  energy  levels,  crack 
branching  is  observed,  Figure  2(b)  (obtained  for  Wo/ Gc  =  13.5). 

While  the  observation  of  consistent  crack  branching  in  successive  simulations  is  an  indication 
of  the  method’s  robustness,  the  detailed  picture  of  energy  dissipation  is  puzzling.  All  simulations 
that  exhibited  branching  showed  an  increase  in  the  amount  of  total  cohesive  energy  dissipated 
as  the  mesh  size  decreased.  This  additional  energy  dissipation  lead  to  decreasing  crack  velocity 
branching  transitions,  although  as  noted  in  Reference  [37]  they  were  consistently  higher  than 
experimentally  observed  values.  Figure  3  shows  energy  dissipation  results  for  the  case  when  the 
ratio  of  strain  energy  stored  in  the  plate  relative  to  that  in  the  cohesive  zone  is  Wo/ Gc  =  13.5. 
A  detailed  examination  of  the  cracked  surfaces  showed  an  increasing  amount  of  microcracking 
and  crack  path  variation  as  the  number  of  nodes  in  the  simulations  increases.  This  trend  has 
been  observed,  albeit  for  a  different  boundary  value  problem,  by  Ruiz  et  al.  [17].  Could  it  be 
that  energy  convergence  may  never  be  attained?  If  numerical  convergence  is  not  an  elusive  goal, 
what  mesh  size  should  be  chosen  so  that  the  amount  of  microcracking  becomes  essentially  mesh 
independent?  All  meshes  but  one  in  Figure  3  resolve  the  cohesize  zone  length  scale.  The  lack 
of  energy  convergence  indicates  that  perhaps  another  criterion  is  needed  for  dynamic  fracture  or 
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fragmentation  simulations.  These  questions  will  be  addressed  in  the  next  section  for  an  idealized 
ID  setting. 


3.  DYNAMIC  FRAGMENTATION  OF  A  CERAMIC  RING 
3.1.  Model  description 

In  order  to  simplify  the  problem,  we  strip  down  the  simulation  capability  and  the  related  complex 
fragmentation  mechanisms  to  the  most  elementary  model.  Our  model  consists  of  an  expanding  ring, 
Figure  4,  which  has  been  discussed  in  earlier  work  [23, 40],  and  was  shown  to  be  equivalent  to  the 
fragmentation  of  a  ID  bar  under  uniform  tension  [41].  As  in  Section  2.4,  the  material  is  assumed 
to  be  linearly  elastic  and  upon  reaching  a  critical  stress,  ID  initially  rigid  cohesive  elements  are 
inserted  to  monitor  multiple  cracks  nucleation  and  interaction.  These  mechanisms  depend  heavily 
on  the  mechanical  properties  of  the  material  and  external  loading.  As  a  crack  grows,  unloading 
stress  waves  are  emitted  and  relieve  the  neighbouring  regions.  The  stress  in  an  unloaded  region 
decreases  with  time  so  that  no  more  cracks  are  initiated,  although  an  existing  crack  may  continue 
to  grow.  If  the  unloading  regions  overlap,  some  nucleated  cracks  may  begin  to  close.  We  assume 
that  stress  waves  pass  through  closed  cracks. 

The  ID  setting  is  convenient  for  two  reasons.  First,  it  allows  the  use  of  the  method  of  charac¬ 
teristics  to  track  stress  wave  interactions  in  the  elastic  material  [40, 42] .  Exact  analytical  solutions 
exist  for  capturing  wave  interactions.  Thus,  the  numerical  part  of  the  method  concentrates  fully 
on  the  non-linear  behaviour  of  cohesive  elements,  which  is  the  focus  of  this  paper.  Nonetheless, 


Broken  cracks 


Closed  cracks; 
Growing  cracks: 

Unloading  Region:  ♦+ 


Figure  4.  Schematic  of  a  ring  fragmenting  under  explosive  loading. 
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the  results  were  also  compared  to  finite-element  calculations  and  no  significant  differences  were 
observed.  Second,  the  restriction  of  the  problem  to  a  ID  setting  permits  the  resolution  of  very  high 
mesh  densities  at  an  affordable  cost.  We  are  then  able  to  assess  unambiguously  if  convergence  can 
be  attained  in  the  context  of  the  cohesive  element  approach. 

Our  model  ceramic  ring  has  a  radius  (Ro)  7.96mm  (the  circumference  of  the  ring,  Lo  =  2nRo ,  is 
50mm).  It  is  made  of  a  material  (similar  to  a  sintered  silicon  nitride)  with  density  2.75  x  103kg/m3, 
elastic  modulus  275  GPa,  fracture  strength  (crc)  300 MPa  and  fracture  energy  (Gc)  lOON/m.  A  1% 
random  variation  around  oc  is  applied  at  all  nodes  to  help  localize  the  process  of  fragmentation. 
The  ring  is  loaded  explosively  with  different  radial  speeds  (vr).  The  equivalent  strain  rate  s  is 
£=vr/Ro.  Each  ring  fragments  through  the  process  described  earlier,  and  fragment  statistics  are 
collected  when  all  the  fragments  are  formed. 

3.2.  Numerical  results 

In  previous  work  we  have  described  the  dependence  of  fragment  size  on  loading  rate  [23,40]. 
Here  we  focus  our  attention  on  the  numerical  convergence  of  the  solution  and  mostly  leave  aside 
physical  discussions. 

Fragmentation  simulations  have  been  conducted  for  13  uniform  mesh  densities  with  num¬ 
ber  of  nodes  ranging  from  128  to  106.  All  these  meshes  have  been  loaded  at  seven  different 
strain  rates  ranging  from  5  x  103  to  5  x  105  s  1 .  Figure  5(a)  represents  the  energy  dissipated  in 
the  process  of  cohesive  fracture.  This  includes  cohesive  elements  that  have  been  fully  damaged 
as  well  as  those  that  are  only  partially  damaged.  For  all  strain  rates  considered  we  observe  a 
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non-monotonic  convergence.  For  clarity,  we  describe  the  results  corresponding  to  the  highest  strain 
rate,  e.g.  5  x  105  s_1.  All  other  strain  rates  follow  a  similar  behaviour,  although  convergence  is 
attained  earlier  and  the  magnitudes  of  energy  dissipation  are  lower  due  to  lesser  amounts  of 
microcracking. 

For  a  strain  rate  of  5  x  105  s-1  (square  symbols  in  Figure  5(a)),  we  see  that  for  coarse  meshes 
the  energy  dissipated  is  proportional  to  the  number  of  nodes  in  the  mesh.  This  trend  is  observed 
roughly  up  to  a  mesh  with  5  x  103  nodes.  This  strong  mesh  dependence  indicates  that  these  mesh 
sizes  under  resolve  the  fragmentation  process.  In  fact,  the  constant  slope  in  Figure  5(a)  implies 
that  all  cohesive  elements  end  up  fully  broken  and  thus  the  fragment  size  depends  directly  on  the 
mesh  size.  Although  Equation  (14)  is  not  directly  applicable  to  a  ID  problem  we  used  it  to  get 
an  estimate  of  the  cohesive  zone  size  for  the  material  properties  under  consideration.  The  vertical 
line  in  Figure  5(a)  represents  the  cohesive  length  scale  based  estimate.  It  indicates  that  meshes 
finer  than  165  nodes  in  the  ID  direction  should  be  able  to  capture  the  cohesive  zone  length  scale. 
Clearly,  for  this  highly  dynamic  problem,  this  mesh  size  estimate  is  not  satisfying,  as  it  cannot 
resolve  the  fragmentation  problem.  In  light  of  this,  there  may  be  a  need  for  alternate  cohesive  zone 
estimates  that  would  include  a  direct  dependence  on  strain  rate.  We  also  note  that  the  dissipated 
energy  increases  up  to  104  nodes  (mesh  A  in  Figure  5).  Although  104  nodes  is  a  relatively  small 
mesh  in  a  ID  setting,  its  2D  equivalent  is  108  nodes  (and  1012  nodes  in  3D).  Many  published 
results  using  the  cohesive  element  approach  in  a  multidimensional  space  have  been  obtained  with 
smaller  meshes.  We  speculate  that  observations  of  energy  increase  with  increasing  number  of 
nodes  in  the  literature  are  due  to  microcracking  being  not  properly  resolved,  although  the  increase 
may  not  be  linearly  proportional  to  the  node  number  as  boundary  conditions  may  deviate  from  the 
simple  uniaxial  state  of  stress  studied  here.  A  relevant  example  is  the  crack  branching  simulations 
in  PMMA  of  this  paper  (Figure  3). 

A  significant  finding  of  this  paper  is  that  for  sufficiently  large  meshes  one  can  properly  resolve 
microcracking.  Indeed,  for  104  nodes  (point  A  in  Figure  5)  a  maximum  in  energy  dissipation  is 
attained  and  beyond  that  point  a  smooth  convergence  is  observed.  For  large  meshes  (106  nodes, 
point  B)  roughly  1000  fragments  are  obtained  (96,  125,  159,  254,  371,  580,  1000,  for  all  seven 
strain  rates  in  increasing  order).  This  finding  is  very  encouraging  as  it  indicates  that  although 
cohesive  elements  are  inserted  at  many  more  nodes  than  necessary  most  are  not  severely  damaged 
and  do  not  impact  significantly  the  energy  balance  of  the  problem.  Despite  this  very  positive 
outcome  for  cohesive  approaches,  it  is  disturbing  to  notice  that  the  same  convergence  may  not 
be  observed  for  meshes  much  smaller  than  106  x  106  =  1012  nodes  in  2D  and  1018  nodes  in  3D, 
which  are  meshes  beyond  any  available  computational  power. 

The  fact  that  convergence  is  non-monotonic  deserves  comments  as  well.  One  may  wonder 
why  a  coarse  mesh  (mesh  A)  leads  to  more  cohesive  cracking  than  a  much  finer  mesh  (mesh 
B).  After  all,  many  more  cohesive  elements  are  present  in  mesh  B  and  these  may  exhibit  more 
microcracking.  Our  initial  attempt  at  explaining  this  trend  consisted  in  checking  if  the  large  time 
step  corresponding  to  the  relatively  coarse  mesh  A  was  too  big  to  unload  sufficiently  slowly  cohesive 
elements  (Equation  (25)).  However,  reducing  the  time  step  by  an  order  of  magnitude  did  not  change 
the  results.  The  problem  is  not  numerical  but  lies  deep  into  the  physics  of  fragmentation  and  more 
precisely  in  the  inherent  randomness  associated  with  generating  fragments.  Figure  6  represents 
fragments  size  distributions  obtained  for  a  ceramic  material  at  various  strain  rates  [41].  We  have 
shown  that  a  simple  universal  law  (Weibull  distribution)  captures  fragment  sizes  for  all  strain 
rates  [42].  The  details  of  this  law  are  beyond  the  scope  of  this  paper  but  the  fact  that  large  and 
small  fragments  exist  at  the  outcome  of  a  fragmentation  event  is  not.  Figures  7(a)  and  (b)  show 
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Figure  6.  PDF  of  fragment  size  for  a  ceramic  ring  under  different  strain  rates  [41]. 


Figure  7.  Fragment  size  distribution  corresponding  to:  (a)  mesh  A;  and  (b)  mesh  B. 


the  fragment  size  distributions  obtained  for  mesh  A  and  mesh  B,  respectively.  The  uniformity  in 
mesh  A  does  not  permit  capturing  the  randomness  of  fragmentation.  Even  though  the  number  of 
nodes  is  relatively  high  most  fragments  vary  little  in  sizes.  The  highly  constrained  fragmentation 
event  ultimately  yields  a  larger  number  of  fragments  (and  therefore  a  higher  dissipated  energy). 
Although  mesh  B  is  still  uniform,  the  number  of  nodes  is  now  sufficiently  large  to  obtain  a  physical 
distribution  of  fragment  sizes,  Figure  7(b). 

Uniform  meshes  severely  constrain  the  fragmentation  event  and  its  desire  to  randomize  fragment 
sizes.  This  explains  the  non-monotonic  convergence  observed  in  Figure  5.  A  natural  question  then 
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Number  of  nodes 


Figure  8.  Effect  of  random  mesh  spacing  on  cohesive  energy  convergence.  The  dashed  line  shows  the 
proposed  fragment- size-based  estimate;  right  of  this  line  convergence  is  attained. 


arises.  Could  random  meshes  lead  to  a  different  convergence  pattern?  The  answer  is  yes  and  this 
in  a  dramatic  way.  The  dotted  lines  in  Figure  8  represent  the  results  obtained  for  random  meshes 
obtained  by  shifting  each  node  by  a  random  amount  in  between  ±0.4/ze,  where  hQ  is  the  mesh 
size  in  the  corresponding  uniform  mesh.  For  this  imposed  degree  of  randomness,  we  obtained  a 
distribution  of  mesh  sizes  with  the  smallest  and  largest  elements  being  around  0.2 hQ  and  1.8 hQ. 
The  convergence  is  now  monotonic  and,  remarkably,  it  is  up  to  two  orders  of  magnitudes  faster. 
For  a  strain  rate  of  5  x  105  s-1,  meshes  with  more  than  104  nodes  seem  to  capture  mostly  the 
fragmentation  event  (Figure  9).  This  relatively  small  number  of  nodes  allows  the  extension  of 
the  analysis  to  higher  dimension  problems.  For  instance,  parallel  fragmentation  simulations  of 
50  mm  x  50  mm  ceramic  plates  under  biaxial  loading,  for  which  108  nodes  may  suffice,  are  within 
sight.  It  is  also  noteworthy  that  the  smoothness  of  the  fragment  size  distribution  obtained  for  a 
mesh  with  106  nodes,  Figure  9,  is  remarkably  improved  compared  to  those  from  uniform  meshes, 
Figure  7(b).  Numerical  modelling  of  fragmentation  has  thus  an  odd  characteristic:  uniform  meshes 
should  be  avoided  at  all  cost. 

Surprisingly,  only  a  slight  degree  of  mesh  randomness  is  necessary  for  improving  numerical 
solution  convergence.  Table  I  lists  the  energy  dissipated  and  number  of  fragments  obtained  by 
imposing  varying  degrees  of  mesh  heterogeneity.  Random  displacements  of  nodes  by  an  amount  in 
between  ±0.05 hQ  suffice  to  improve  convergence  by  roughly  two  orders  of  magnitude.  It  is  however 
possible  that  an  optimum  mesh  size  distribution  may  be  found  to  achieve  fastest  convergence.  There, 
fragmentation  analytical  studies  based  on  the  principle  of  maximum  entropy  [43-45]  may  prove 
useful  to  determine  the  mesh  density  functional. 


Copyright  ©  2006  John  Wiley  &  Sons,  Ltd. 


Int.  J.  Numer.  Meth.  Engng  2007;  69:484-503 
DOI:  10.1002/nme 


Frequtincr  aium 


THE  COHESIVE  ELEMENT  APPROACH  TO  DYNAMIC  FRAGMENTATION 


499 


Figure  9.  Mesh  dependence  of  fragment  size  distribution  (random 
mesh  spacing)  for  5  x  105  s-1  strain  rate. 


Table  I.  Effect  of  different  perturbation  factors  on  energy  dissipated  and 
number  of  fragments  (for  5  x  105  s-1). 


Perturbation  factor 
(as  %  of  hG) 

Energy  dissipated 
(mJ) 

#  fragments 

0 

393.4 

2723 

5 

152.1 

1093 

10 

147.2 

1065 

15 

143.8 

1042 

20 

142.0 

1047 

25 

139.7 

1038 

30 

136.2 

1015 

35 

136.4 

1018 

40 

134.7 

994 

Copyright  ©  2006  John  Wiley  &  Sons,  Ltd. 


Int.  J.  Numer.  Meth.  Engng  2007;  69:484-503 
DOI:  10.1002/nme 


500 


J.  F.  MOLINARI  ET  AL. 


As  discussed  earlier,  at  least  in  a  ID  setup,  cohesive  zone  estimates  such  as  Equation  (14),  lead 
to  meshes  too  coarse  to  capture  all  microcracking  events.  It  is  therefore  of  interest  to  establish  an 
alternate  estimate,  even  if  empirical.  A  key  Ending  of  our  previous  work  [42]  was  that  a  universal 
law  may  be  used  to  capture  the  dependence  of  the  average  fragment  size  with  strain  rate.  This 
estimate  was  proposed  as  an  alternate  to  Grady’s  energy-based  model  [46],  which  does  not  take 
into  account  dynamical  effects.  The  proposed  equation,  which  captures  accurately  a  wide  range  of 
ceramic  materials,  is 


(27) 


where  s  is  the  average  fragment  size,  s  is  the  strain  rate,  £o  =  (oc/E)/t$  is  the  characteristic  strain 
rate  (Equation  (18)  [36]),  and  so  =  cto  is  the  characteristic  fragment  size.  These  latter  parameters 
both  depend  on  the  characteristic  time  to  =  pc5c/2oc  (Equation  (13)). 

Many  fragments  end  up  being  smaller  than  the  average  fragment  size.  In  a  ID  setting,  a  mesh 
size  that  is  roughly  one  order  of  magnitude  smaller  than  the  average  size  leads  to  converged  results. 
Our  mesh  size  estimate,  expanded  to  explicitly  present  all  the  material  parameters,  is 


he=0.25(EGc/a2c) 


1+4.5 


£ 

co\!  E2Gc 


(28) 


The  number  of  nodes  corresponding  to  Equation  (28)  is  shown  by  the  dashed  line  in  Figure  8. 
It  seems  to  be  a  reasonable  indicator  of  the  occurrence  of  convergence  for  the  range  of  strain 
rates  studied.  Other  estimates,  perhaps  incorporating  time  scales  such  as  Equation  (25),  may  be 
proposed.  Based  on  these  estimates,  converged  2D  fragmentation  simulations  are  within  reach. 
3D  simulations,  on  the  other  hand,  remain  perhaps  too  computationally  intensive  to  achieve  fully 
converged  solutions.  There,  cohesive  approaches  would  be  well  served  to  take  into  account  the 
poor  resolution  of  microcracking  events.  Rate-dependent  cohesive  laws,  for  which  Gc  becomes 
an  increasing  function  of  strain  rate,  are  promising  in  this  regard  [37,47,48].  In  addition,  a 
worthy  research  direction  consists  in  checking,  for  dynamic  fragmentation  problems,  convergence 
properties  of  other  numerical  approaches.  In  this  regard  partition  of  unity  methods  [49]  as  well  as 
discrete  and  meshfree  methods  [50, 51],  may  lead  to  faster  convergence  as  crack  locations  are  not 
constrained  to  appear  at  element  boundaries.  Initially,  these  cross  comparisons  could  be  conducted 
in  a  ID  setting. 


4.  CONCLUSIONS 

The  main  theme  of  this  paper  has  been  energy  convergence  of  the  cohesive  element  approach 
for  the  fragmentation  analysis  of  a  linear  elastic  material.  While  the  discussion  has  not  addressed 
crack  path  convergence,  substantial  computational  evidence  demonstrates,  for  the  first  time,  that 
the  cohesive  element  method  converges  in  an  energetic  sense,  at  least  in  a  simple  one-dimensional 
setting.  Microcracking  events  and  the  ensuing  fragment  sizes  distributions  are  statistically  mesh 
independent  for  sufficiently  fine  meshes.  Remarkably,  convergence  was  attained  up  to  two  orders 
of  magnitude  earlier  for  random  meshes  than  for  uniform  meshes,  and  that  even  for  very  small 
random  perturbations.  It  should  be  emphasized  that  standard  cohesive  zone  size  estimates  under 
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resolve  the  number  of  nodes  necessary  for  attaining  mesh  independence  in  dynamic  fragmentation 
simulations.  This  may  explain  why  energy  convergence  of  the  cohesive  element  approach  had 
not  been  previously  observed.  We  have  proposed  a  simple  mesh  density  estimate  based  on  the 
dependence  of  the  average  fragment  size  on  the  strain  rate  and  elastic  and  fracture  material 
properties.  This  estimate  provides  a  clear  roadmap  for  extending  the  findings  of  this  paper  to  more 
complex  loading  conditions  including  biaxial  loading  of  ceramic  plates. 
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